Load libraries

library(tidyverse)
library(scran)
library(scater)
library(BiocParallel)
library(gridExtra)
library(ComplexHeatmap)
library(circlize)
library(ggbeeswarm)
library(readxl)
library(parallel)
library(edgeR)

ncores <- 10
mcparam <- MulticoreParam(workers=ncores)
register(mcparam)

Define variables

opt <- list()
opt$sce <- "data/submission/TFlow_split_complete.RDS"
opt$types <- "data/submission/TFlow_celltypes.csv"
opt$plot <- "plots/Fig6/"
opt$drugscreen <- "data/submission/drugScreens_pseudo.RDS"

## Set theme
lgd <-  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), 
        text = element_text(size = 15),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent',colour = NA), 
        strip.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

Load data

## Load single-cell experiment
sce <- readRDS(opt$sce)
sce.unstim <- sce[["unstim"]]
sce.stim <- sce[["stim"]]
rm(sce) # free-up the memory

## Load cell type annotation
types <- read.csv(opt$types) %>%
  mutate(uniqueBarcode = as.character(uniqueBarcode)) %>%
  mutate(celltype_broad = ifelse(celltype %in% c("CD16- NK-cell", "CD16+ NK-cell"), "NK-Cells", celltype)) %>%
  mutate(celltype_broad = ifelse(celltype %in% c("T-reg", "gd T-cell", "naive CD4 T-cell", "naive CD8 T-cell", "CD8 T-emra", "CD8 effector memory T-cell", 
                                                 "CD8 central memory T-cell", "CD7+ memory CD4 T-emra", "CD7+ effector memory CD4 T-cell", 
                                                 "CD7+ central memory CD4 T-cell", "CD7- memory CD4 T-emra", "CD7- effector memory CD4 T-cell", 
                                                 "CD7- central memory CD4 T-cell"), "T-Cells", celltype_broad)) %>%
  mutate(celltype_broad = ifelse(celltype %in% c("naive B-cell", "memory B-cell"), "B-Cells", celltype_broad)) %>%
  mutate(celltype = ifelse(celltype %in% c("CD7- effector memory CD4 T-cell", "CD7+ effector memory CD4 T-cell"), "CD4 EM T-cell", celltype), 
         celltype = ifelse(celltype %in% c("CD7- central memory CD4 T-cell", "CD7+ central memory CD4 T-cell"), "CD4 CM T-cell", celltype), 
         celltype = ifelse(celltype %in% c("CD8 central memory T-cell"), "CD8 CM T-cell", celltype),
         celltype = ifelse(celltype %in% c("CD8 effector memory T-cell"), "CD8 EM T-cell", celltype),
         celltype = ifelse(celltype %in% c("CD7- effector memory CD4 T-cell", "CD7+ effector memory CD4 T-cell"), "CD4 EM T-cell", celltype),
         celltype = ifelse(celltype %in% c("CD7- memory CD4 T-emra", "CD7+ memory CD4 T-emra"), "CD4 T-emra", celltype))

## Drug screen data
drugList <- readRDS(opt$drugscreen)

## Define drug combinations
drugComb <- c("Birinapant|Necrostatin-1 (25)", "Birinapant|QVD-Oph (25)", "Birinapant|Necrostatin-1 (25)|QVD-Oph (25)", 
              "Birinapant|Necrostatin-1 (12.5)", "Birinapant|QVD-Oph (12.5)", "Birinapant|Necrostatin-1 (12.5)|QVD-Oph (12.5)",
              "GDC-0152|Necrostatin-1 (25)", "GDC-0152|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (25)|QVD-Oph (25)", "GDC-0152|Necrostatin-1 (12.5)", 
              "GDC-0152|QVD-Oph (12.5)", "GDC-0152|Necrostatin-1 (12.5)|QVD-Oph (12.5)", "Birinapant|Ipatasertib (2)", 
              "Birinapant|Ruxolitinib (2)", "Birinapant|Dacinostat (0.04)", "Birinapant|Venetoclax (0.04)",
              "Birinapant|Bafilomycin_A1 (0.08)", "Birinapant|NSA (0.8)", "Birinapant|NSA (2)", 
              "Birinapant|NSA (5)", "Birinapant|NSA (12.5)", "Birinapant|NSA (0.8)|QVD-Oph (25)", 
              "Birinapant|NSA (2)|QVD-Oph (25)", "Birinapant|NSA (5)|QVD-Oph (25)", "Birinapant|NSA (12.5)|QVD-Oph (25)", 
              "Birinapant + Ipatasertib 2µM", "Birinpant + Ruxolitinib 2µM", "Birinapant + Bafilomycin A1 0.08µM", 
              "Birinapant + NSA 0.8µM",
              "Birinapant + NSA 0.8µM + QVD-Oph 25µM", "Birinapant + NSA 2µM",  "Birinapant + NSA 5µM", 
              "Birinapant + NSA 2µM + QVD-Oph 25µM", "Birinapant + NSA 5μM + QVD-Oph 25μM", 
              "DMSO", "empty", "Necrostatin-1|QVD-Oph")

Analysis

Overview - Stimulated

ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45
il2 <- 1.5
tnf  <- 1.5
ifng <- 1.25
ki67 <- 1.5
gzmb <- 1.5
gmcsf <- 1.5
il10 <- 1.5
gzmb <- 1.5

plotTab <- ggcells(sce.stim, features = c("CD3", "FSCA", "SSCA", "Ki67", "Apotracker", "LD", "cleaved_caspase_3", 
                                     "IL2", "IL10", "TNF", "IFNg", "GMCSF"), exprs_values = "exprs")[[1]]

## Randomly sample to reduce overplotting
set.seed(100)
plotTab <- plotTab[sample(1:nrow(plotTab), nrow(plotTab)), ]

## Plot Treatment
p <- plotTab %>%
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Treatment)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd + #scale_colour_viridis_c() +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("DMSO_stim" = "grey70", "Birinapant_low_stim" = "#AD1457"
                               ))
p

ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_treatment_stim.png"), 
       height = 15, width = 15, units = "cm")

## Plot Diagnosis
p <- plotTab %>%
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Diagnosis_simple)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd + #scale_colour_viridis_c() +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6","healthy" = "#FFA000"
                               ))
p

ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_diagnosis_stim.png"), 
       height = 15, width = 15, units = "cm")

## Plot celltype
p <- plotTab %>%
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
  #mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd + #scale_colour_viridis_c() +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6", "T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
                                     "B-Cells" = "#FFA000"#, "monocyte" = "#66BB6A""#F44336"
                               ))
p

ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_celltype_stim.png"), 
       height = 15, width = 15, units = "cm")

p1 <- plotTab %>% 
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype)) %>%
  mutate(celltype_broad = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
         celltype_broad = ifelse(FSCA < 750000, "dead", celltype_broad),
         celltype_broad = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad),
         celltype_broad = ifelse(celltype == "dead" & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad)) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("live" = "#A5D6A7","apoptotic" = "#FFA000", 
                                "dead" = "#F44336"
                               ))
p1

ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celldeath_stim.png"), 
       height = 15, width = 15, units = "cm")

Overview - Unstimulated

plotTab <- ggcells(sce.unstim, features = c("CD3", "LD", "Apotracker", "cleaved_caspase_3", "FSCA", "CD25"), exprs_values = "exprs")[[1]]

## Randomly sample to reduce overplotting
set.seed(100)
plotTab <- plotTab[sample(1:nrow(plotTab), nrow(plotTab)), ]

ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45

## Plot cell types
p1 <- plotTab %>% 
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype), 
         Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6",
                                "T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
                                "B-Cells" = "#FFA000"
                               ))
p1

ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celltypes.png"), 
       height = 15, width = 15, units = "cm")

## Plot cell death
p1 <- plotTab %>% 
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype), 
         Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
  mutate(celltype_broad = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"), 
         celltype_broad = ifelse(FSCA < 750000, "dead", celltype_broad),
         celltype_broad = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad),
         celltype_broad = ifelse(celltype == "dead" & cleaved_caspase_3 > cl3, "apoptotic", celltype_broad)) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = celltype_broad)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("live" = "#A5D6A7","apoptotic" = "#FFA000", 
                                "dead" = "#F44336"
                               ))
p1

ggsave(plot = p1, filename = paste0(opt$plot, "TFlow_umap_celldeath.png"), 
       height = 15, width = 15, units = "cm")


## Plot treatments
p <- plotTab %>%
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype), 
         Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Treatment)) +
  geom_point(size = 0.05) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("Birinapant" = "#AD1457", "QVDOph" = "#64B5F6", 
                                 "Birinapant|QVDOph" = "#43A047", "Selinexor" = "#FFD45F"
                               ))
p

ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_treatments.png"), 
       height = 15, width = 15, units = "cm")

## Plot Diagnosis
p <- plotTab %>%
  filter(!is.na(UMAP_FSCA.1)) %>%
  left_join(types, by = "uniqueBarcode") %>%
  filter(!celltype %in% c("debris"), !is.na(celltype), 
         Treatment %in% c("DMSO", "Birinapant", "Selinexor", "QVDOph", "Birinapant|QVDOph")) %>%
  ggplot(aes(x = UMAP_FSCA.1, y = UMAP_FSCA.2, colour = Diagnosis_simple)) +
  geom_point(size = 0.1) +
  xlab("UMAP1") + ylab("UMAP2") +
  lgd + #scale_colour_viridis_c() +
  theme_void() +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        legend.position = "none") +
  scale_colour_manual(values = c("T-LGL" = "#A5D6A7", "T-PLL" = "#64B5F6","healthy" = "#FFA000"
                               ))
p

ggsave(plot = p, filename = paste0(opt$plot, "TFlow_umap_diagnosis.png"), 
       height = 15, width = 15, units = "cm")

Show cell death frequency relative to DMSO

plotTab <- ggcells(sce.unstim, features = c("CD3", "LD", "Apotracker", "cleaved_caspase_3", "FSCA", "CD25", "CD127"), exprs_values = "exprs")[[1]]

# Rare combinations of celltype x patient x treatment will be exclduded from the analysis, as this makes it more susceptible to noise. 
countSub <- plotTab %>% 
  left_join(types, by = "uniqueBarcode") %>%
  dplyr::count(patientID, Treatment, celltype_broad) %>%
  filter(n < 100) %>%
  mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
  pull(patTreatType) %>% unique()

plotTab %>%
  select(patientID, Treatment) %>% unique() %>%
  dplyr::count(Treatment)
##           Treatment  n
## 1        Birinapant 15
## 2 Birinapant|QVDOph 15
## 3              DMSO 15
## 4            QVDOph 15
## 5         Selinexor 15
compDrug <- c("Selinexor", "Birinapant_low", "Birinapant", 
              "Birinapant|QVDOph", "QVDOph",
              "DMSO_stim", "Birinapant_low_stim")

ld <- 1.5
cl3 <- 2.5
cl3small <- 2.25
apo <- 1.45

testTab <- mclapply(unique(plotTab$patientID), mc.cores = ncores, function(x) {
  testTab <- plotTab %>% filter(patientID == x) %>%
  left_join(types, by = "uniqueBarcode") %>% 
  mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
  mutate(cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"), 
                  cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
         cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
         cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
  mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
  filter(!patTreatType %in% countSub) %>% # filter only for combinations with enough cells
  group_by(patientID, celltype_broad, Treatment) %>%
  mutate(nCells = length(uniqueBarcode)) %>%
  ungroup() %>%
  group_by(patientID, celltype_broad, cellDeath, Treatment) %>%
  mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
  select(patientID, Diagnosis_simple, celltype_broad, Treatment, cellDeath, rel) %>% unique() %>%
  pivot_wider(names_from = "cellDeath", values_from = "rel") %>%
  mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
  pivot_longer(cols = c("live", "apoptotic", "dead"), 
               names_to = "cellDeath", values_to = "rel") %>%
  mutate(rel = ifelse(is.na(rel), 0, rel))

}) %>% bind_rows()

## Normalize to DMSO
viabTab <- mclapply(unique(testTab$patientID), mc.cores = ncores, function(x) {
  lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
    lapply(c("live", "apoptotic", "dead"), function(y) {
      message(x)
      message(z)
      message(y)
      
      testTab <- testTab %>% filter(patientID == x, cellDeath == y, celltype_broad == z)
      
      dmso.val <- unique(testTab[testTab$Treatment == "DMSO", ]$rel)
      
      if(length(dmso.val == 1)) {
        testTab %>% mutate(viabNorm = rel / dmso.val, 
                           viabDiff = rel - dmso.val)
      } else {
        testTab %>% mutate(viabNorm = NA)
      }
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows()

viabTab.back <- viabTab

compDrug <- c("Selinexor", "Birinapant_low", "Birinapant", 
              "Birinapant|Q-VD-OPh", "Q-VD-OPh",
              "DMSO_stim", "Birinapant_low_stim")

# Out of range observations
viabTab <- viabTab %>% mutate(viabNorm = ifelse(cellDeath == "live" & viabNorm > 1.25, 1.25, viabNorm), 
                              viabNorm = ifelse(cellDeath == "apoptotic" & viabNorm > 5, 5, viabNorm), 
                              viabNorm = ifelse(cellDeath == "dead" & viabNorm > 5, 5, viabNorm)) %>%
  mutate(Treatment = str_replace(string = Treatment, pattern = "QVDOph", replacement = "Q-VD-OPh"))

viabTabOOR <- viabTab[viabTab$cellDeath == "live" & viabTab$viabNorm == 1.25 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]

## Plot normalized viability
p <- viabTab %>%
  filter(cellDeath == "live", !is.na(viabNorm)) %>%
  filter(Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh"
                          )) %>%
  mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
  ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
  geom_boxplot(alpha = 0.6) +
  geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 1.25 & viabTab$cellDeath == "live" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
  geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +
  #geom_text(aes(label = patientID)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  xlab("") + ylab("Living Cells") +
  lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), 
              legend.position = "none", 
              text = element_text(size = 20), 
              panel.border = element_rect(linewidth = 1, fill = "transparent")) +
  facet_wrap(. ~ celltype_broad, nrow = 1) +
  scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6", 
                               "Birinapant|Q-VD-OPh" = "#43A047"))
p

ggsave(plot = p, filename = paste0(opt$plot, "live_viabNorm_celltypes.png"), 
       height = 14, width = 30, units = "cm")

## Show the cell death type
viabTabOOR <- viabTab[viabTab$cellDeath == "apoptotic" & viabTab$viabNorm == 5 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]


## Apoptotic
p <- viabTab %>%
  filter(cellDeath == "apoptotic",
         !is.na(viabNorm)) %>%
  filter(!Treatment %in% c("DMSO")) %>%
  mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
  ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
  geom_boxplot(alpha = 0.6) +
  geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 5 & viabTab$cellDeath == "apoptotic" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
  geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +

  geom_hline(yintercept = 1, linetype = "dashed") +
  xlab("") + ylab("Apoptotic Cells") +
  lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), 
              legend.position = "none", 
              text = element_text(size = 20), 
              panel.border = element_rect(linewidth = 1, fill = "transparent")) +
  facet_wrap(. ~ celltype_broad, nrow = 1) +
  scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6", 
                               "Birinapant|Q-VD-OPh" = "#43A047"))
p

ggsave(plot = p, filename = paste0(opt$plot, "apoptotic_viabNorm_celltypes.png"), 
       height = 14, width = 30, units = "cm")

## Non-apoptotic
viabTabOOR <- viabTab[viabTab$cellDeath == "dead" & viabTab$viabNorm == 5 & viabTab$Treatment %in% c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"), ]
viabTabOOR <- viabTabOOR[!is.na(viabTabOOR$patientID), ]

p <- viabTab %>%
  filter(cellDeath == "dead", #celltype_broad %in% c("lymphoma"), 
         !is.na(viabNorm)) %>%
  filter(!Treatment %in% c("DMSO")) %>%
  mutate(Treatment = factor(Treatment, levels = compDrug)) %>%
  ggplot(aes(x = Treatment, y = viabNorm, fill = Treatment)) +
  geom_boxplot(alpha = 0.6) +
  geom_beeswarm(cex = 3.25, size = 2.5, shape = 21, data = viabTab[viabTab$viabNorm < 5 & viabTab$cellDeath == "dead" & viabTab$Treatment %in% c("Selinexor", "Birinapant", "Q-VD-OPh", "Birinapant|Q-VD-OPh") & !is.na(viabTab$viabNorm), ]) +
  geom_beeswarm(cex = 3, size = 2.25, shape = 24, data = viabTabOOR) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  xlab("") + ylab("Non-Apoptotic Dying Cells") +
  lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), 
              legend.position = "none", 
              text = element_text(size = 20), 
              panel.border = element_rect(linewidth = 1, fill = "transparent")) +
  facet_wrap(. ~ celltype_broad, nrow = 1) +
  scale_fill_manual(values = c("Birinapant" = "#AD1457", "Selinexor" = "#FFD45F", "Q-VD-OPh" = "#64B5F6", 
                               "Birinapant|Q-VD-OPh" = "#43A047"))
p

ggsave(plot = p, filename = paste0(opt$plot, "non_apoptotic_viabNorm_celltypes.png"), 
       height = 14, width = 30, units = "cm")

Perform T-test

de.res <- lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
    lapply(c("live", "apoptotic", "dead"), function(y) {
      lapply(c("Birinapant", "Selinexor", "Q-VD-OPh", "Birinapant|Q-VD-OPh"
               ), function(x) {
        message(z)
        message(y)
        message(x)
        testTab <- viabTab %>% filter(celltype_broad == z, cellDeath == y, Treatment %in% c("DMSO", x)) %>%
          mutate(Treatment = factor(Treatment, levels = c("DMSO", x))) %>%
          select(-viabNorm, -viabDiff) %>%
          pivot_wider(names_from = "Treatment", values_from = "rel")

        if(nrow(testTab) >= 1) {
          colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
          #testTab <- testTab[complete.cases(testTab), ]

          res <- t.test(testTab$compTreat, testTab$DMSO, alternative = "two.sided",
                 var.equal = TRUE, paired = TRUE)

          resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
        }

        # t.test(testTab$compTreat, alternative = "less",
        #        var.equal = TRUE, mu = 1)

      }) %>% bind_rows()
    }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"), 
                              p.adj = round(p.adj, 3), 
                              diff = round(diff, 2))
rownames(de.res) <- NULL
de.res
##    celltype_broad cellDeath            p  diff           Treatment p.adj
## 1  Lymphoma Cells      live 2.127039e-03 -0.10          Birinapant 0.009
## 2  Lymphoma Cells      live 1.180840e-02 -0.23           Selinexor 0.023
## 3  Lymphoma Cells      live 1.074801e-01  0.06            Q-VD-OPh 0.129
## 4  Lymphoma Cells      live 3.270611e-03 -0.44 Birinapant|Q-VD-OPh 0.010
## 5  Lymphoma Cells apoptotic 7.680785e-01  0.01          Birinapant 0.801
## 6  Lymphoma Cells apoptotic 2.414977e-02  0.12           Selinexor 0.040
## 7  Lymphoma Cells apoptotic 1.243085e-02 -0.12            Q-VD-OPh 0.023
## 8  Lymphoma Cells apoptotic 8.204419e-01  0.01 Birinapant|Q-VD-OPh 0.838
## 9  Lymphoma Cells      dead 3.120977e-04  0.09          Birinapant 0.003
## 10 Lymphoma Cells      dead 1.288734e-02  0.11           Selinexor 0.023
## 11 Lymphoma Cells      dead 1.035210e-01  0.06            Q-VD-OPh 0.127
## 12 Lymphoma Cells      dead 5.487479e-04  0.43 Birinapant|Q-VD-OPh 0.004
## 13        B-Cells      live 1.138709e-02 -0.14          Birinapant 0.023
## 14        B-Cells      live 1.194643e-03 -0.35           Selinexor 0.006
## 15        B-Cells      live 2.115535e-02  0.04            Q-VD-OPh 0.036
## 16        B-Cells      live 2.706574e-02 -0.20 Birinapant|Q-VD-OPh 0.042
## 17        B-Cells apoptotic 6.132582e-02  0.05          Birinapant 0.089
## 18        B-Cells apoptotic 4.441851e-03  0.30           Selinexor 0.013
## 19        B-Cells apoptotic 9.382079e-04 -0.17            Q-VD-OPh 0.006
## 20        B-Cells apoptotic 1.934377e-03 -0.13 Birinapant|Q-VD-OPh 0.008
## 21        B-Cells      dead 4.671281e-02  0.09          Birinapant 0.070
## 22        B-Cells      dead 9.250892e-02  0.05           Selinexor 0.123
## 23        B-Cells      dead 9.290184e-03  0.13            Q-VD-OPh 0.019
## 24        B-Cells      dead 6.522731e-03  0.33 Birinapant|Q-VD-OPh 0.016
## 25        T-Cells      live 9.522270e-02 -0.08          Birinapant 0.124
## 26        T-Cells      live 9.186048e-04 -0.22           Selinexor 0.006
## 27        T-Cells      live 1.880645e-03  0.10            Q-VD-OPh 0.008
## 28        T-Cells      live 3.326217e-01 -0.09 Birinapant|Q-VD-OPh 0.355
## 29        T-Cells apoptotic 1.007910e-01  0.04          Birinapant 0.127
## 30        T-Cells apoptotic 2.171900e-04  0.17           Selinexor 0.003
## 31        T-Cells apoptotic 7.737700e-05 -0.15            Q-VD-OPh 0.002
## 32        T-Cells apoptotic 7.070642e-03 -0.10 Birinapant|Q-VD-OPh 0.016
## 33        T-Cells      dead 3.108288e-01  0.04          Birinapant 0.339
## 34        T-Cells      dead 1.552538e-01  0.05           Selinexor 0.177
## 35        T-Cells      dead 7.977198e-02  0.05            Q-VD-OPh 0.109
## 36        T-Cells      dead 7.603422e-02  0.18 Birinapant|Q-VD-OPh 0.107
## 37       NK-Cells      live 1.326006e-04 -0.08          Birinapant 0.002
## 38       NK-Cells      live 5.257889e-03 -0.14           Selinexor 0.013
## 39       NK-Cells      live 5.290138e-03  0.10            Q-VD-OPh 0.013
## 40       NK-Cells      live 1.125779e-01 -0.07 Birinapant|Q-VD-OPh 0.132
## 41       NK-Cells apoptotic 3.057807e-03  0.06          Birinapant 0.010
## 42       NK-Cells apoptotic 1.221021e-05  0.10           Selinexor 0.001
## 43       NK-Cells apoptotic 2.714987e-03 -0.15            Q-VD-OPh 0.010
## 44       NK-Cells apoptotic 9.244809e-01  0.00 Birinapant|Q-VD-OPh 0.924
## 45       NK-Cells      dead 2.703958e-02  0.02          Birinapant 0.042
## 46       NK-Cells      dead 2.209660e-01  0.04           Selinexor 0.247
## 47       NK-Cells      dead 5.101816e-03  0.06            Q-VD-OPh 0.013
## 48       NK-Cells      dead 8.600919e-03  0.08 Birinapant|Q-VD-OPh 0.019
## Combination index
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
  lapply(c("live")[[1]], function(x) {
  lapply(c("Birinapant|QVDOph"), function(y) {
    message(z)
    message(x)
    message(y)
    testTab <- viabTab.back %>% 
      filter(celltype_broad == z, cellDeath == x, Treatment %in% c("DMSO", "QVDOph", "Birinapant", y)) %>%
      select(-rel, -viabDiff) %>% 
      pivot_wider(names_from = "Treatment", values_from = "viabNorm")
            
      colnames(testTab)[[which(colnames(testTab) == y)]] <- "compTreat"

      if(y == "Birinapant|QVDOph"){
          testTab <- testTab %>% mutate(expViab = Birinapant * QVDOph)
        }
      ## Remove samples with two NAs
      testTab <- testTab %>% mutate(#compTreat = ifelse(is.na(compTreat), 0, compTreat),
                                    CI = compTreat / expViab) %>% filter(!is.na(CI))

      res <- t.test(testTab$expViab, testTab$CI, alternative = "two.sided", var.equal = TRUE, paired = TRUE)
      resTab <- data.frame(celltype_broad = z, cellDeath = x, Treatment = y, diff = res$estimate, p = res$p.value,
                           meanCI = mean(testTab$CI))

      }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"),
                              p.adj = round(p.adj, 3), 
                              meanCI = round(meanCI, 2))
##                     celltype_broad cellDeath         Treatment         diff
## mean difference...1 Lymphoma Cells      live Birinapant|QVDOph  0.447124843
## mean difference...2        B-Cells      live Birinapant|QVDOph -0.576507228
## mean difference...3        T-Cells      live Birinapant|QVDOph  0.052968832
## mean difference...4       NK-Cells      live Birinapant|QVDOph -0.008682567
##                              p meanCI p.adj
## mean difference...1 0.01631067   0.39 0.065
## mean difference...2 0.41684653   1.38 0.834
## mean difference...3 0.68971695   0.88 0.898
## mean difference...4 0.89829792   0.95 0.898

Show the proportions among live cells

countSub.sub <- plotTab %>% 
  left_join(types, by = "uniqueBarcode") %>%
  mutate(celltype = ifelse(celltype %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype)) %>%
  mutate( cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
                  cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
         cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
         cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
  filter(cellDeath == "live") %>%
  dplyr::count(patientID, Treatment, celltype) %>%
  filter(n < 100) %>%
  mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype)) %>%
  pull(patTreatType) %>% unique()

#### Now look at the proportions of healthy cells
testTab <- mclapply(unique(plotTab[plotTab$Diagnosis_simple == "healthy", ]$patientID), mc.cores = ncores, function(x) {
  plotTab %>% filter(patientID == x) %>%
  left_join(types, by = "uniqueBarcode") %>% 
  filter(!celltype_broad %in% c("debris")) %>%
  mutate(celltype = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype), 
         celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
  mutate(cellDeath = ifelse(LD > ld | cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
                  cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
         cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
         cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
  filter(cellDeath == "live") %>%
  mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype)) %>%
  filter(!patTreatType %in% countSub.sub) %>% # filter only for combinations with enough cells
  group_by(patientID, Treatment) %>%
  mutate(nCells = length(uniqueBarcode)) %>%
  ungroup() %>%
  group_by(patientID, celltype, Treatment) %>%
  mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
  select(patientID, Diagnosis_simple, celltype_broad, celltype, Treatment, nCellsSub, nCells, cellDeath, rel) %>% unique() 

}) %>% bind_rows()

freqTab <- mclapply(unique(testTab$patientID), mc.cores = ncores, function(x) {
  lapply(unique(testTab$celltype), function(z) {
    lapply(c("live"), function(y) {
      message(x)
      message(z)
      message(y)
      
      testTab <- testTab %>% filter(patientID == x, cellDeath == y, celltype == z)
      
      dmso.val <- unique(testTab[testTab$Treatment == "DMSO", ]$rel)
      
      if(length(dmso.val == 1)) {
        testTab %>% mutate(viabNorm = rel / dmso.val, 
                           viabDiff = rel - dmso.val)
      } else {
        testTab %>% mutate(viabNorm = NA)
      }
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows()

lapply(c("Birinapant", "Selinexor"), function(x) {
  lapply(c("T-Cells"), function(y) {
    df <- freqTab %>%
      filter(cellDeath == "live", !celltype_broad %in% c("debris"),
             Treatment == x, viabNorm != Inf) %>%
      group_by(Treatment, celltype) %>%
      summarise(meanViab = mean(viabNorm, na.rm = TRUE), 
                sdViab = sd(viabNorm, na.rm = TRUE) / sqrt(length(viabNorm))) %>%
      ungroup() 
    
    cellOrder <- df %>%
      arrange(desc(meanViab)) %>%
      pull(celltype) %>%
      unique()
    
    p <- df %>%
      mutate(celltype = factor(celltype, levels = cellOrder),
             celltype_broad = ifelse(celltype %in% c("CD16- NK-cell", "CD16+ NK-cell"), "NK-Cells", NA),
             celltype_broad = ifelse(celltype %in% c("naive B-cell", "memory B-cell"), "B-Cells", celltype_broad),
             celltype_broad = ifelse(celltype %in% c("T-PLL", "T-LGL", "Lymphoma Cells"), "Lymphoma Cells", celltype_broad),
             celltype_broad = ifelse(celltype %in% c("CD8 EM T-cell", "CD8 CM T-cell",
                                                     "CD8 T-emra", "naive CD8 T-cell", "gd T-cell", "CD4 CM T-cell",
                                                     "naive CD4 T-cell", "CD4 EM T-cell", "T-reg",
                                                     "CD4 T-emra"), "T-Cells", celltype_broad)) %>%
     ggplot(aes(x = celltype, y = meanViab, fill = celltype_broad)) +
        geom_bar(stat = "identity", colour = "black", alpha = 0.6) +
        geom_beeswarm(aes(x = celltype, y = viabNorm), size = 3, shape = 21, cex = 2, data = freqTab[freqTab$Treatment == x, ]) +
        geom_errorbar(aes(ymin = meanViab, ymax = meanViab + sdViab), width = 0.5) +
        geom_hline(yintercept = 1, linetype = "dashed") +
        xlab("") +
      ylab("Frequency Rel. to DMSO") + ggtitle(paste0(x, " Effect on Healthy Donors (Living Cells)")) +
      lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
                  legend.position = "none",
                  panel.border = element_rect(linewidth = 1, fill = "transparent"),
                  axis.text = element_text(size = 15),
                  text = element_text(size = 15)) +
      scale_fill_manual(values = c("Lymphoma Cells" = "#64B5F6", "T-Cells" = "#0D47A1", "NK-Cells" = "#AB47BC",
                                    "B-Cells" = "#FFA000"
                               ))
    
    ggsave(plot = p, filename = paste0(opt$plot, x, "_diff_effect_all_subtypes_bar.png"),
           height = 14, width = 21, units = "cm")
    p
    
   })
})
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
## [[1]]
## [[1]][[1]]

## 
## 
## [[2]]
## [[2]][[1]]

Compute the fraction of proliferating cells

plotTab <- ggcells(sce.stim, features = c("FSCA", "Ki67", "Apotracker", "LD", "cleaved_caspase_3", "CD45RA", "CCR7",
                                     "IL2", "IL10", "TNF", "IFNg", "GMCSF", "GranzymeB", "CD16", "CD56", "CD3", "CD19", "CD4", "CD8", "CD27", "TCRgd"), exprs_values = "exprs")[[1]]

## Get live cells
brc.live <- plotTab %>%
    left_join(types, by = "uniqueBarcode") %>% 
   mutate(cellDeath = ifelse(LD > ld |  cleaved_caspase_3 > cl3 | Apotracker > apo, "dead", "live"),
         cellDeath = ifelse(FSCA < 750000, "dead", cellDeath),
         cellDeath = ifelse(FSCA < 1500000 & cleaved_caspase_3 > cl3small, "apoptotic", cellDeath),
         cellDeath = ifelse(cellDeath == "dead" & cleaved_caspase_3 > cl3, "apoptotic", cellDeath)) %>%
  filter(cellDeath %in% c("live")) %>%
  pull(uniqueBarcode) %>% unique()

## Get only samples with decent cell counts
countSub <- plotTab %>% 
  filter(uniqueBarcode %in% brc.live) %>%
  left_join(types, by = "uniqueBarcode") %>%
  dplyr::count(patientID, Treatment, celltype_broad) %>%
  filter(n < 100) %>%
  mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
  pull(patTreatType) %>% unique()

compDrug <- c("Selinexor", "Birinapant_low", "Birinapant", "Birinapant|QVDOph", "QVDOph",
              "DMSO_stim", "Birinapant_low_stim")

il2 <- 1.5
tnf  <- 1.5
ifng <- 1.25
ki67 <- 1.5
gzmb <- 1.5
gmcsf <- 1.5
il10 <- 1.5

testTab <- mclapply(c("IL2", "Ki67", "TNF", "IFNg", "GMCSF", "IL10"), mc.cores = ncores, 
                  function(x) {
  if(x == "IL2") {cut.off <- il2}
  if(x == "Ki67") {cut.off <- ki67}
  if(x == "TNF") {cut.off <- tnf}
  if(x == "IFNg") {cut.off <- ifng}
  if(x == "GMCSF") {cut.off <- gmcsf}
  if(x == "IL10") {cut.off <- il10}
  
  plotTab %>%
    left_join(types, by = "uniqueBarcode") %>% 
    mutate(celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
    filter(FSCA < 3000000, uniqueBarcode %in% brc.live) %>%
    pivot_longer(cols = x, names_to = "var", values_to = "expr") %>%
    mutate(cellCytokine = ifelse(expr > cut.off, "pos", "neg")) %>%
    mutate(patTreatType = paste0(patientID, "_", Treatment, "_", celltype_broad)) %>%
    filter(!patTreatType %in% countSub) %>% # filter only for combinations with enough cells
    group_by(patientID, celltype_broad, Treatment) %>%
    mutate(nCells = length(uniqueBarcode)) %>%
    ungroup() %>%
    group_by(patientID, celltype_broad, cellCytokine, Treatment) %>%
    mutate(nCellsSub = length(uniqueBarcode), rel = nCellsSub / nCells) %>%
    select(patientID, Diagnosis_simple, celltype_broad, Treatment, cellCytokine, rel) %>% unique() %>%
    pivot_wider(names_from = "cellCytokine", values_from = "rel") %>%
    mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
    pivot_longer(cols = c("pos", "neg"), 
                 names_to = "cellCytokine", values_to = "rel") %>%
    mutate(rel = ifelse(is.na(rel), 0, rel)) %>%
    mutate(cytokine = x)
}) %>% bind_rows()

## Normalize to DMSO
viabTab <- mclapply(c("IL2", "Ki67", "TNF", "IFNg", "GMCSF", "IL10"), mc.cores = ncores, function(y) {
  lapply(unique(testTab$patientID), function(x) {
  lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {

      message(x)
      message(z)
      #message(y)
      
      testTab <- testTab %>% filter(patientID == x, cytokine == y, cellCytokine == "pos", celltype_broad == z)
      
      dmso.val <- unique(testTab[testTab$Treatment == "DMSO_stim", ]$rel)
      
      if(length(dmso.val == 1)) {
        testTab %>% mutate(viabNorm = rel / dmso.val)
      } else {
        testTab %>% mutate(viabNorm = NA)
      }
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows()


## Are there differences between stimulated and unstimulated DMSO
lapply(c("Ki67"), function(x) {
  
  p <- viabTab %>% filter(cytokine == x) %>%
  filter(Treatment %in% c("DMSO_stim", "Birinapant_low_stim")) %>%
  mutate(Treatment = as.character(Treatment), 
         Treatment = ifelse(Treatment == "DMSO_stim", "DMSO", Treatment)) %>%
  mutate(Treatment = ifelse(Treatment == "Birinapant_low_stim", "Birinapant", Treatment)) %>%
  mutate(Treatment = factor(Treatment, levels = c("DMSO", compDrug))) %>%
  ggplot(aes(x = Treatment, y = rel, fill = Treatment)) +
  geom_line(linewidth = 0.5, aes(group = patientID)) +
  geom_boxplot(alpha = 0.6) +
  geom_beeswarm(cex = 3, size = 3, shape = 21) +
  xlab("") + ylab(paste0("Frac. of ", x, " Expressing Cells")) +
  ggtitle("Birinapant Effect on Proliferation") +
  lgd + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1), 
              legend.position = "none", 
              text = element_text(size = 17.5), 
              panel.border = element_rect(linewidth = 1, fill = "transparent"), 
              axis.text = element_text(size = 17.5)) +
  scale_fill_manual(values = c("Birinapant" = "#AD1457", "DMSO" = "grey70")) +
  facet_wrap(. ~ celltype_broad, nrow = 1)
  
  if(x == "Ki67") {
    ggsave(plot = p, filename = paste0(opt$plot, "TFlow_celltype_prol.png"),
        height = 12, width = 30, units = "cm" 
       )
  }
p

})
## [[1]]

Perform paired t-test

lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
    lapply(c("IL2", "IL10", "TNF", "IFNg", "GMCSF", "Ki67")[[6]], function(y) {
      lapply(c("Birinapant_low_stim"), function(x) {
        message(z)
        message(y)
        message(x)
        testTab <- viabTab %>% filter(celltype_broad == z, cellCytokine == "pos", cytokine == y, Treatment %in% c("DMSO_stim", x)) %>%
          mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", x))) %>%
          select(-viabNorm) %>%
          pivot_wider(names_from = "Treatment", values_from = "rel")

        if(nrow(testTab) >= 1) {
          colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
          #testTab <- testTab[complete.cases(testTab), ]

          res <- t.test(testTab$compTreat, testTab$DMSO_stim, alternative = "two.sided",
                 var.equal = TRUE, paired = TRUE)

          resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
        }

        # t.test(testTab$compTreat, alternative = "less",
        #        var.equal = TRUE, mu = 1)
        
      }) %>% bind_rows()
    }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"), 
                              p.adj = round(p.adj, 3), 
                              diff = round(diff, 2))
##                     celltype_broad cellDeath          p  diff
## mean difference...1 Lymphoma Cells      Ki67 0.01894206 -0.02
## mean difference...2        B-Cells      Ki67 0.91953751  0.00
## mean difference...3        T-Cells      Ki67 0.09410205 -0.01
## mean difference...4       NK-Cells      Ki67 0.52465156  0.00
##                               Treatment p.adj
## mean difference...1 Birinapant_low_stim 0.076
## mean difference...2 Birinapant_low_stim 0.920
## mean difference...3 Birinapant_low_stim 0.188
## mean difference...4 Birinapant_low_stim 0.700
lapply(c("Lymphoma Cells", "B-Cells", "T-Cells", "NK-Cells"), function(z) {
    lapply(c("IL2", "TNF", "IFNg", "GMCSF"), function(y) {
      lapply(c("Birinapant_low_stim"), function(x) {
        message(z)
        message(y)
        message(x)
        testTab <- viabTab %>% filter(celltype_broad == z, cellCytokine == "pos", cytokine == y, Treatment %in% c("DMSO_stim", x)) %>%
          mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", x))) %>%
          select(-viabNorm) %>%
          pivot_wider(names_from = "Treatment", values_from = "rel")

        if(nrow(testTab) >= 1) {
          colnames(testTab)[[which(colnames(testTab) == x)]] <- "compTreat"
          #testTab <- testTab[complete.cases(testTab), ]

          res <- t.test(testTab$compTreat, testTab$DMSO_stim, alternative = "two.sided",
                 var.equal = TRUE, paired = TRUE)

          resTab <- data.frame(celltype_broad = z, cellDeath = y, p = res$p.value, diff = res$estimate, Treatment = x)
        }

        # t.test(testTab$compTreat, alternative = "less",
        #        var.equal = TRUE, mu = 1)
        
      }) %>% bind_rows()
    }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"), 
                              p.adj = round(p.adj, 3), 
                              diff = round(diff, 2))
##                      celltype_broad cellDeath           p  diff
## mean difference...1  Lymphoma Cells       IL2 0.340359218  0.01
## mean difference...2  Lymphoma Cells       TNF 0.057969995 -0.03
## mean difference...3  Lymphoma Cells      IFNg 0.350525204 -0.02
## mean difference...4  Lymphoma Cells     GMCSF 0.022189361  0.01
## mean difference...5         B-Cells       IL2 0.651648886  0.00
## mean difference...6         B-Cells       TNF 0.004911755 -0.05
## mean difference...7         B-Cells      IFNg 0.373888420  0.00
## mean difference...8         B-Cells     GMCSF 0.912209574  0.00
## mean difference...9         T-Cells       IL2 0.170823607  0.06
## mean difference...10        T-Cells       TNF 0.252360329 -0.03
## mean difference...11        T-Cells      IFNg 0.423812004 -0.01
## mean difference...12        T-Cells     GMCSF 0.027517255  0.01
## mean difference...13       NK-Cells       IL2 0.362919164  0.00
## mean difference...14       NK-Cells       TNF 0.019890190  0.12
## mean difference...15       NK-Cells      IFNg 0.103392260  0.13
## mean difference...16       NK-Cells     GMCSF 0.016077270  0.17
##                                Treatment p.adj
## mean difference...1  Birinapant_low_stim 0.460
## mean difference...2  Birinapant_low_stim 0.155
## mean difference...3  Birinapant_low_stim 0.460
## mean difference...4  Birinapant_low_stim 0.088
## mean difference...5  Birinapant_low_stim 0.695
## mean difference...6  Birinapant_low_stim 0.079
## mean difference...7  Birinapant_low_stim 0.460
## mean difference...8  Birinapant_low_stim 0.912
## mean difference...9  Birinapant_low_stim 0.342
## mean difference...10 Birinapant_low_stim 0.449
## mean difference...11 Birinapant_low_stim 0.484
## mean difference...12 Birinapant_low_stim 0.088
## mean difference...13 Birinapant_low_stim 0.460
## mean difference...14 Birinapant_low_stim 0.088
## mean difference...15 Birinapant_low_stim 0.236
## mean difference...16 Birinapant_low_stim 0.088

Test for differential expression

sce.x <- sce.stim
compDrug <- c("Birinapant_low_stim")


colData(sce.x) <- colData(sce.x) %>%
  data.frame() %>%
  left_join(types, by = "uniqueBarcode") %>%
  mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", compDrug)), 
         cluster_id = celltype_broad, 
         celltype_broad = ifelse(celltype_broad %in% c("T-PLL", "T-LGL"), "Lymphoma Cells", celltype_broad)) %>%
  DataFrame()

testTreat <- c("Birinapant_low_stim")
cellSub <- c("Lymphoma Cells", "T-Cells", "B-Cells", "NK-Cells")
adtSub <- c("IL2", "GMCSF", "TNF", "IFNg")

## For some reason it still thinks there are NAs
sce.x <- sce.x[, !is.na(sce.x$celltype_broad)]
sce.x <- sce.x[, sce.x$uniqueBarcode %in% brc.live]

de.res <- lapply(testTreat, function(z) {
  message(paste0("Fitting model for ", z))
  lapply(cellSub, function(x) {
    message(paste0("Celltype ", x))
    
    sce.x <- sce.x[, sce.x$celltype_broad == x & sce.x$uniqueBarcode %in% brc.live]
    
    ## Set-up model formulas
    meta.df <- unique(dplyr::select(data.frame(colData(sce.x)), patientID, Treatment, Stimulation))
    design <- stats::model.matrix(~ patientID + Treatment, meta.df)
    frml <- diffcyt::createFormula(meta.df, cols_fixed = c("patientID", "Treatment"))
    
    ## Get median counts
    summed <- summarizeAssayByGroup(sce.x, ids = colData(sce.x)[,c("patientID", "Treatment")],
                                    statistics = "median", assay.type = "exprs") # get raw counts and sum them up
    colData(summed)$Treatment <- factor(colData(summed)$Treatment, levels = c("DMSO_stim", compDrug))
    
    mat <- assay(summed, "median") %>% as.matrix()
    colnames(mat) <- paste0(summed$patientID, ".", summed$Treatment)
    
    ## For every marker assemble the corresponding data frame and fit a linear model
    lapply(adtSub, function(u) {
      df <- data.frame(medianExpr = mat[u, ], patTreat = colnames(mat), ncells = summed$ncells) %>%
        separate(col = "patTreat", into = c("patientID", "Treatment"), sep = "\\.")
      
      testTab <- frml$data %>% left_join(df, by = c("patientID", "Treatment")) %>%
        dplyr::rename(y = medianExpr) %>%
        mutate(Treatment = factor(Treatment, levels = c("DMSO_stim", compDrug))) %>%
        filter(!ncells < 100)
      testTab
      
      ## Fit the linear model
      fit <- lm(frml$formula, data = testTab, weights = testTab$ncells)

      ## Fit contrasts
      cntr <- rep(0, times = length(names(fit$coefficients)))
      cntr[[which(names(fit$coefficients) == paste0("Treatment", z))]] <- 1
      cntr <- diffcyt::createContrast(cntr) %>% t()

      res <- multcomp::glht(fit, cntr)

      ## Assemble output
      resTab <- data.frame(celltype = x, Treatment = z, adt = u, estimate = summary(res)$test$coefficient, p = summary(res)$test$pvalues)
    }) %>% bind_rows()
  }) %>% bind_rows()
}) %>% bind_rows() %>% mutate(p.adj = p.adjust(p, method = "BH"))
de.res %>%
  mutate(estimate = round(estimate, 2), 
         p.adj = round(p.adj, 3))
##              celltype           Treatment   adt estimate           p p.adj
## 1...1  Lymphoma Cells Birinapant_low_stim   IL2     0.10 0.278423532 0.318
## 1...2  Lymphoma Cells Birinapant_low_stim GMCSF     0.06 0.080219214 0.227
## 1...3  Lymphoma Cells Birinapant_low_stim   TNF    -0.10 0.130927389 0.246
## 1...4  Lymphoma Cells Birinapant_low_stim  IFNg    -0.04 0.318886545 0.340
## 1...5         T-Cells Birinapant_low_stim   IL2     0.32 0.085254230 0.227
## 1...6         T-Cells Birinapant_low_stim GMCSF    -0.03 0.258601059 0.318
## 1...7         T-Cells Birinapant_low_stim   TNF    -0.11 0.252182432 0.318
## 1...8         T-Cells Birinapant_low_stim  IFNg    -0.07 0.179300332 0.287
## 1...9         B-Cells Birinapant_low_stim   IL2     0.00 0.983271082 0.983
## 1...10        B-Cells Birinapant_low_stim GMCSF     0.02 0.138553496 0.246
## 1...11        B-Cells Birinapant_low_stim   TNF    -0.12 0.004702416 0.073
## 1...12        B-Cells Birinapant_low_stim  IFNg     0.02 0.039789897 0.159
## 1...13       NK-Cells Birinapant_low_stim   IL2    -0.02 0.251544631 0.318
## 1...14       NK-Cells Birinapant_low_stim GMCSF     0.64 0.009137436 0.073
## 1...15       NK-Cells Birinapant_low_stim   TNF     0.37 0.017227989 0.092
## 1...16       NK-Cells Birinapant_low_stim  IFNg     0.23 0.129153616 0.246
p <- de.res %>%
  mutate(pSign = -log10(p.adj) * sign(estimate), 
         label = ifelse(p.adj < 0.1, "*", "")) %>%
  ggplot(aes(x = celltype, y = adt)) +
  geom_tile(aes(fill = pSign), colour = "black") +
  geom_text(aes(label = label), size = 6.5) +
  xlab("") + ylab("") + ggtitle("Birinapant Effect on Cytokine\nExpression (< 10% FDR)") +
  lgd + theme(panel.border = element_rect(colour = "black", fill=NA, linewidth = 1), 
              text = element_text(size = 17.5),
              axis.text = element_text(size = 17.5),
              axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  scale_fill_gradient2(low = "#1976D2", high = "#F44336",  
                       name = "-Log10(p.adj) \nwith Direction")
p

ggsave(plot = p, filename = paste0(opt$plot, "cytokine_heatmap.png"), 
       height = 15, width = 17.5, units = "cm")

Output session info

sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] parallel  grid      stats4    stats     graphics  grDevices datasets 
##  [8] utils     methods   base     
## 
## other attached packages:
##  [1] edgeR_4.0.16                limma_3.58.1               
##  [3] readxl_1.4.3                ggbeeswarm_0.7.2           
##  [5] circlize_0.4.16             ComplexHeatmap_2.18.0      
##  [7] gridExtra_2.3               BiocParallel_1.36.0        
##  [9] scater_1.30.1               scran_1.30.2               
## [11] scuttle_1.12.0              SingleCellExperiment_1.24.0
## [13] SummarizedExperiment_1.32.0 Biobase_2.62.0             
## [15] GenomicRanges_1.54.1        GenomeInfoDb_1.38.8        
## [17] IRanges_2.36.0              S4Vectors_0.40.2           
## [19] BiocGenerics_0.48.1         MatrixGenerics_1.14.0      
## [21] matrixStats_1.3.0           lubridate_1.9.3            
## [23] forcats_1.0.0               stringr_1.5.1              
## [25] dplyr_1.1.4                 purrr_1.0.2                
## [27] readr_2.1.5                 tidyr_1.3.1                
## [29] tibble_3.2.1                ggplot2_3.5.1              
## [31] tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2               bitops_1.0-7               
##   [3] cellranger_1.1.0            polyclip_1.10-7            
##   [5] XML_3.99-0.17               lifecycle_1.0.4            
##   [7] rstatix_0.7.2               doParallel_1.0.17          
##   [9] lattice_0.21-9              MASS_7.3-60                
##  [11] backports_1.5.0             magrittr_2.0.3             
##  [13] sass_0.4.9                  rmarkdown_2.27             
##  [15] jquerylib_0.1.4             yaml_2.3.9                 
##  [17] metapod_1.10.1              minqa_1.2.7                
##  [19] RColorBrewer_1.1-3          ConsensusClusterPlus_1.66.0
##  [21] multcomp_1.4-25             abind_1.4-5                
##  [23] zlibbioc_1.48.2             Rtsne_0.17                 
##  [25] RCurl_1.98-1.14             TH.data_1.1-2              
##  [27] tweenr_2.0.3                sandwich_3.1-0             
##  [29] GenomeInfoDbData_1.2.11     ggrepel_0.9.5              
##  [31] irlba_2.3.5.1               dqrng_0.4.1                
##  [33] DelayedMatrixStats_1.24.0   codetools_0.2-19           
##  [35] DelayedArray_0.28.0         ggforce_0.4.2              
##  [37] tidyselect_1.2.1            shape_1.4.6.1              
##  [39] farver_2.1.2                lme4_1.1-35.5              
##  [41] ScaledMatrix_1.10.0         viridis_0.6.5              
##  [43] jsonlite_1.8.8              GetoptLong_1.0.5           
##  [45] BiocNeighbors_1.20.2        survival_3.5-7             
##  [47] iterators_1.0.14            systemfonts_1.1.0          
##  [49] foreach_1.5.2               tools_4.3.2                
##  [51] ggnewscale_0.5.0            ragg_1.3.2                 
##  [53] Rcpp_1.0.12                 glue_1.7.0                 
##  [55] SparseArray_1.2.4           xfun_0.45                  
##  [57] withr_3.0.0                 fastmap_1.2.0              
##  [59] boot_1.3-28.1               bluster_1.12.0             
##  [61] fansi_1.0.6                 digest_0.6.36              
##  [63] rsvd_1.0.5                  timechange_0.3.0           
##  [65] R6_2.5.1                    textshaping_0.4.0          
##  [67] colorspace_2.1-0            utf8_1.2.4                 
##  [69] generics_0.1.3              renv_1.0.7                 
##  [71] S4Arrays_1.2.1              pkgconfig_2.0.3            
##  [73] gtable_0.3.5                RProtoBufLib_2.14.1        
##  [75] XVector_0.42.0              diffcyt_1.22.1             
##  [77] htmltools_0.5.8.1           carData_3.0-5              
##  [79] clue_0.3-65                 scales_1.3.0               
##  [81] png_0.1-8                   colorRamps_2.3.4           
##  [83] knitr_1.48                  rstudioapi_0.16.0          
##  [85] reshape2_1.4.4              tzdb_0.4.0                 
##  [87] rjson_0.2.21                nlme_3.1-163               
##  [89] nloptr_2.1.1                cachem_1.1.0               
##  [91] zoo_1.8-12                  GlobalOptions_0.1.2        
##  [93] vipor_0.4.7                 pillar_1.9.0               
##  [95] vctrs_0.6.5                 ggpubr_0.6.0               
##  [97] BiocSingular_1.18.0         car_3.1-2                  
##  [99] cytolib_2.14.1              beachmat_2.18.1            
## [101] cluster_2.1.4               beeswarm_0.4.0             
## [103] evaluate_0.24.0             mvtnorm_1.2-5              
## [105] cli_3.6.3                   locfit_1.5-9.10            
## [107] compiler_4.3.2              rlang_1.1.4                
## [109] crayon_1.5.3                ggsignif_0.6.4             
## [111] labeling_0.4.3              FlowSOM_2.10.0             
## [113] plyr_1.8.9                  flowCore_2.14.2            
## [115] stringi_1.8.4               viridisLite_0.4.2          
## [117] munsell_0.5.1               Matrix_1.6-1.1             
## [119] hms_1.1.3                   sparseMatrixStats_1.14.0   
## [121] statmod_1.5.0               highr_0.11                 
## [123] igraph_2.0.3                broom_1.0.6                
## [125] bslib_0.7.0